Last updated: 2023-08-11

Checks: 7 0

Knit directory: scSeq_Hefendehl/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20220131)

The command set.seed(20220131) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 9d29eff

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 9d29eff. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated: 


Ignored files:
    Ignored:    .RData
    Ignored:    .Rproj.user/
    Ignored:    data/20230203/
    Ignored:    data/Alignments/
    Ignored:    data/Cluster.markers_02-02-2023.xlsx
    Ignored:    data/DAM_genelists.RData
    Ignored:    data/DEG_control_vs_stroke.xlsx
    Ignored:    data/ReloadAllData_Hefendehl_Stroke_Dec'21.RData
    Ignored:    data/Sample_Tables/
    Ignored:    data/Stroke-SS2_SCT-harmony_Hefendehl.html
    Ignored:    data/counts.csv
    Ignored:    data/genecounts.csv
    Ignored:    data/microglia-SCT_02-02-23.rds
    Ignored:    data/microglia_protein.rds
    Ignored:    data/samples.integrated.RData
    Ignored:    data/seu-SCT-harmony_02-02-23.rds
    Ignored:    data/seu_filtered_nonorm_02-02-23.rds
    Ignored:    data/tx2genes.csv
    Ignored:    output/CrossFCplot.pdf
    Ignored:    output/CrossFCplot.svg
    Ignored:    output/integrated_analysis.RData
    Ignored:    output/integrated_analysis.RDataTmp

Untracked files:
    Untracked:  shinySecret.R

Unstaged changes:
    Modified:   geneviewer/Dataset.RData
    Modified:   geneviewer/rsconnect/shinyapps.io/molgenlab/geneviewer.dcf
    Modified:   wflowhelper.R

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/01_Biostat.Rmd) and HTML (docs/01_Biostat.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd 9d29eff Andreas Chiocchetti 2023-08-11 workflowr::wflow_publish(c("./docs/", "./analysis/", "./code/*",
html 9d29eff Andreas Chiocchetti 2023-08-11 workflowr::wflow_publish(c("./docs/", "./analysis/", "./code/*",
html 87fbacb Andreas Chiocchetti 2023-08-07 Build site.
Rmd e103c55 Andreas Chiocchetti 2023-08-07 workflowr::wflow_publish(c("./docs/", "./analysis/", "./code/*",
Rmd d2da683 Andreas Chiocchetti 2023-08-06 update markdowns
html d2da683 Andreas Chiocchetti 2023-08-06 update markdowns
html 6d3e973 Andreas Chiocchetti 2023-07-14 Build site.
Rmd d5c8125 Andreas Chiocchetti 2023-07-14 workflowr::wflow_publish(c("./docs/", "./analysis/", "./code/*",
html d5c8125 Andreas Chiocchetti 2023-07-14 workflowr::wflow_publish(c("./docs/", "./analysis/", "./code/*",
html 982b9f6 Andreas Chiocchetti 2023-02-20 Build site.
Rmd 438f7ad Andreas Chiocchetti 2023-02-20 workflowr::wflow_publish(c("./docs/", "./analysis/", "./code/*",
html 0e20f22 Andreas Chiocchetti 2023-02-16 Build site.
Rmd eaa373d Andreas Chiocchetti 2023-02-16 workflowr::wflow_publish(c("./docs/", "./analysis/", "./code/*",
html eaa373d Andreas Chiocchetti 2023-02-16 workflowr::wflow_publish(c("./docs/", "./analysis/", "./code/*",
html 7764484 achiocch 2022-06-27 Build site.
html 14e1ee0 achiocch 2022-06-27 Build site.
Rmd 5bba014 achiocch 2022-06-27 wflow_publish(c("analysis/", "docs/", "code/*"))
html 5bba014 achiocch 2022-06-27 wflow_publish(c("analysis/", "docs/", "code/*"))
html 3451e15 achiocch 2022-06-22 Build site.
Rmd cfbfcf6 achiocch 2022-06-22 wflow_publish(c("analysis/", "docs/", "code/*"))
html cfbfcf6 achiocch 2022-06-22 wflow_publish(c("analysis/", "docs/", "code/*"))
html faf20e9 achiocch 2022-05-16 Build site.
Rmd 8791434 achiocch 2022-05-16 wflow_publish(c("analysis/", "docs/", "code/*"))
html 8791434 achiocch 2022-05-16 wflow_publish(c("analysis/", "docs/", "code/*"))
html 51b35c8 achiocch 2022-04-26 Build site.
Rmd 3277c08 achiocch 2022-04-26 wflow_publish(c("analysis/", "docs/", "code/*"))
html 41d6cd0 achiocch 2022-04-26 Build site.
Rmd 7d20262 achiocch 2022-04-26 wflow_publish(c("analysis/", "docs/", "code/*"))
Rmd 8e44cfc achiocch 2022-04-22 wflow_publish(c("analysis/", "docs/", "code/*"))
html 8e44cfc achiocch 2022-04-22 wflow_publish(c("analysis/", "docs/", "code/*"))
html 9b778b1 achiocch 2022-04-19 Build site.
Rmd d2ed458 achiocch 2022-04-19 wflow_publish(c("analysis/", "docs/", "code/*"))
html d2ed458 achiocch 2022-04-19 wflow_publish(c("analysis/", "docs/", "code/*"))
html 55b3c82 achiocch 2022-04-08 wflow_publish(c("analysis/", "docs/", "code/*"))
html cf395f4 achiocch 2022-03-30 Build site.
Rmd 55e7738 achiocch 2022-03-30 wflow_publish(c("analysis/", "docs/", "code/*"))
html 55e7738 achiocch 2022-03-30 wflow_publish(c("analysis/", "docs/", "code/*"))
html c998828 achiocch 2022-03-30 Build site.
Rmd f21fdf7 achiocch 2022-03-30 wflow_publish(c("analysis/", "docs/", "code/*"))
html f21fdf7 achiocch 2022-03-30 wflow_publish(c("analysis/", "docs/", "code/*"))
html 87439a3 achiocch 2022-02-21 Build site.
Rmd ef56d05 achiocch 2022-02-21 wflow_publish(c("analysis/", "docs/", "code/*"))
html d6f2105 achiocch 2022-02-21 Build site.
Rmd b832d1c achiocch 2022-02-21 wflow_publish(c("analysis/", "docs/", "code/*"))
Rmd 95fb700 achiocch 2022-02-21 wflow_publish(c("analysis/", "docs/", "code/*"))
html 95fb700 achiocch 2022-02-21 wflow_publish(c("analysis/", "docs/", "code/*"))
html 79dd0b1 achiocch 2022-02-21 Build site.
Rmd 73aabec achiocch 2022-02-21 wflow_publish(c("analysis/", "docs/", "code/*"))
html ded601d achiocch 2022-02-03 Build site.
Rmd 858d646 achiocch 2022-02-03 wflow_publish(c("analysis/", "docs/", "code/*"))
html 858d646 achiocch 2022-02-03 wflow_publish(c("analysis/", "docs/", "code/*"))
html 50fe211 achiocch 2022-02-03 Build site.
Rmd 64cd37c achiocch 2022-02-03 wflow_publish(c("analysis/", "docs/", "code/*"))
html 5fbdba4 achiocch 2022-02-02 Build site.
Rmd 636a117 achiocch 2022-02-02 wflow_publish(c("analysis/", "docs/", "code/*"))
html 636a117 achiocch 2022-02-02 wflow_publish(c("analysis/", "docs/", "code/*"))
Rmd 9044798 Andreas Geburtig-Chiocchetti 2022-02-02 workflowr automated
Rmd 57c9e70 achiocch 2022-02-02 adds initial html files
html 57c9e70 achiocch 2022-02-02 adds initial html files
Rmd ddb6373 achiocch 2022-02-02 first commit 
#get sample data
samples.integrated@meta.data %>% as.data.frame() -> samplemeta


spacial_data= samples.microglia@meta.data %>% as.data.frame() %>% select(contains("spatial"))

samplemeta <- cbind(samplemeta, spacial_data[rownames(samplemeta),])

# convert to correct data type 

# define genotype as is
samplemeta$Genotype_corr = factor(samplemeta$Genotype=="WT", levels=c(F,T), labels = c("APPPS1+", "WT"))
samplemeta$Genotype_corr = relevel(samplemeta$Genotype_corr, ref="WT")

samplemeta$methoxy = factor(samplemeta$Methoxy=="pos", levels=c(T,F), labels = c("MX04+", "MX04-"))

samplemeta$Treatment = as.factor(samplemeta$Treatment)
samplemeta$Treatment = relevel(samplemeta$Treatment, ref="Ctrl")
samplemeta$Mouse_ID = as.factor(samplemeta$Mouse_ID)
samplemeta$Sex = as.factor(samplemeta$Gender)
samplemeta$Brain_region = as.factor(samplemeta$Brain_region)
samplemeta$Celltype = as.factor(samplemeta$Celltype)

samples.integrated@meta.data = samplemeta

nCells=nrow(samplemeta)
nMice=nlevels(samplemeta$Mouse_ID)
nCelltypes=nlevels(samplemeta$Celltype)

Sample descriptive:

Data contains a total of 562 Cells from 9.

Cells per Mouse

Cells per Strain

table(Mouse_ID=samplemeta$Mouse_ID) %>% as.data.frame() %>% display_tab()

Cells per Genotype

table(Genotype=samplemeta$Genotype_corr, Treatment=samplemeta$Treatment) %>% as.data.frame() %>% display_tab()
table(Genotype=samplemeta$Genotype_corr, Celltype=samplemeta$Celltype) %>% 
  as.data.frame() %>% display_tab()
## add spatial cell types 
table(Genotype=samplemeta$Genotype_corr, Sp_Controls=samplemeta$Control_spatial_binary) %>% 
  as.data.frame() %>% display_tab()
table(Genotype=samplemeta$Genotype_corr, Sp_Stroke=samplemeta$Stroke_spatial_binary) %>% 
  as.data.frame() %>% display_tab()
table(Genotype=samplemeta$Genotype_corr, 
      Celltype=samplemeta$Celltype, Treatment=samplemeta$Treatment
) %>% 
  as.data.frame() %>% display_tab()
samplemeta$condition=paste0(samplemeta$Genotype_corr,"_", samplemeta$Treatment)
samples.integrated@meta.data = samplemeta


variables=c("Celltype","Sex", "Age","Genotype_corr", "Treatment", "condition","Phase","Control_spatial_binary","Stroke_spatial_binary", "Brain_region","nCount_RNA","pseudoaligned_reads", "percent.mt", "percent.rp", "Mouse_ID")

Descriptive stats across Cell type

res = compareGroups(Celltype~., data = samplemeta[,variables], max.ylev = 10)
#summary(res)
export_table <- createTable(res)
options(width = 10000)
export2md(export_table)
Summary descriptives table by groups of `Celltype’
Microglia_1 Microglia_2 Microglia_3 Microglia_4 Granulocytes Microglia_5 Mono/Mac p.overall
N=372 N=86 N=57 N=8 N=4 N=15 N=20
Sex: 0.276
f 73 (19.6%) 20 (23.3%) 17 (29.8%) 2 (25.0%) 2 (50.0%) 5 (33.3%) 5 (25.0%)
m 299 (80.4%) 66 (76.7%) 40 (70.2%) 6 (75.0%) 2 (50.0%) 10 (66.7%) 15 (75.0%)
Age: .
37 weeks 226 (60.8%) 56 (65.1%) 41 (71.9%) 6 (75.0%) 2 (50.0%) 7 (46.7%) 11 (55.0%)
38 weeks 44 (11.8%) 8 (9.30%) 13 (22.8%) 2 (25.0%) 2 (50.0%) 4 (26.7%) 3 (15.0%)
40 weeks 102 (27.4%) 22 (25.6%) 3 (5.26%) 0 (0.00%) 0 (0.00%) 4 (26.7%) 6 (30.0%)
Genotype_corr: .
WT 170 (45.7%) 43 (50.0%) 42 (73.7%) 8 (100%) 4 (100%) 8 (53.3%) 15 (75.0%)
APPPS1+ 202 (54.3%) 43 (50.0%) 15 (26.3%) 0 (0.00%) 0 (0.00%) 7 (46.7%) 5 (25.0%)
Treatment: .
Ctrl 272 (73.1%) 56 (65.1%) 52 (91.2%) 8 (100%) 4 (100%) 11 (73.3%) 8 (40.0%)
Stroke 100 (26.9%) 30 (34.9%) 5 (8.77%) 0 (0.00%) 0 (0.00%) 4 (26.7%) 12 (60.0%)
condition: .
APPPS1+_Ctrl 164 (44.1%) 28 (32.6%) 10 (17.5%) 0 (0.00%) 0 (0.00%) 6 (40.0%) 1 (5.00%)
APPPS1+_Stroke 38 (10.2%) 15 (17.4%) 5 (8.77%) 0 (0.00%) 0 (0.00%) 1 (6.67%) 4 (20.0%)
WT_Ctrl 108 (29.0%) 28 (32.6%) 42 (73.7%) 8 (100%) 4 (100%) 5 (33.3%) 7 (35.0%)
WT_Stroke 62 (16.7%) 15 (17.4%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 3 (20.0%) 8 (40.0%)
Phase: .
G1 159 (42.7%) 36 (41.9%) 19 (33.3%) 0 (0.00%) 2 (50.0%) 10 (66.7%) 9 (45.0%)
G2M 95 (25.5%) 18 (20.9%) 17 (29.8%) 2 (25.0%) 2 (50.0%) 2 (13.3%) 6 (30.0%)
S 118 (31.7%) 32 (37.2%) 21 (36.8%) 6 (75.0%) 0 (0.00%) 3 (20.0%) 5 (25.0%)
Control_spatial_binary: 0.080
FALSE 338 (91.1%) 75 (87.2%) 57 (100%) 8 (100%) 4 (100%) 15 (100%) 8 (100%)
TRUE 33 (8.89%) 11 (12.8%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 0 (0.00%)
Stroke_spatial_binary: .
FALSE 240 (64.7%) 55 (64.0%) 53 (93.0%) 8 (100%) 4 (100%) 11 (73.3%) 8 (100%)
TRUE 131 (35.3%) 31 (36.0%) 4 (7.02%) 0 (0.00%) 0 (0.00%) 4 (26.7%) 0 (0.00%)
Brain_region: .
Cortex 272 (73.1%) 56 (65.1%) 52 (91.2%) 8 (100%) 4 (100%) 11 (73.3%) 8 (40.0%)
Lesion 100 (26.9%) 30 (34.9%) 5 (8.77%) 0 (0.00%) 0 (0.00%) 4 (26.7%) 12 (60.0%)
nCount_RNA 158416 (103720) 194411 (106729) 164921 (154435) 217767 (205467) 204215 (86217) 165673 (126013) 162375 (90426) 0.154
pseudoaligned_reads 158947 (104087) 194836 (106825) 165764 (154492) 218138 (205598) 204425 (86225) 166053 (126105) 163005 (90604) 0.160
percent.mt 1.90 (1.35) 2.14 (1.78) 1.44 (1.07) 1.66 (1.04) 0.98 (1.13) 2.01 (1.01) 2.04 (1.27) 0.089
percent.rp 2.21 (1.29) 3.05 (2.17) 1.48 (1.34) 0.98 (0.41) 1.02 (0.50) 2.03 (1.25) 3.16 (2.30) <0.001
Mouse_ID: .
256#1022 33 (8.87%) 9 (10.5%) 13 (22.8%) 2 (25.0%) 0 (0.00%) 1 (6.67%) 2 (10.0%)
256#1023 31 (8.33%) 11 (12.8%) 16 (28.1%) 4 (50.0%) 2 (50.0%) 0 (0.00%) 2 (10.0%)
364#469 44 (11.8%) 8 (9.30%) 13 (22.8%) 2 (25.0%) 2 (50.0%) 4 (26.7%) 3 (15.0%)
386 16 (4.30%) 4 (4.65%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 4 (20.0%)
387 9 (2.42%) 3 (3.49%) 1 (1.75%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 2 (10.0%)
388 77 (20.7%) 15 (17.4%) 2 (3.51%) 0 (0.00%) 0 (0.00%) 4 (26.7%) 0 (0.00%)
409 46 (12.4%) 11 (12.8%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 3 (20.0%) 4 (20.0%)
457 87 (23.4%) 13 (15.1%) 8 (14.0%) 0 (0.00%) 0 (0.00%) 2 (13.3%) 1 (5.00%)
461 29 (7.80%) 12 (14.0%) 4 (7.02%) 0 (0.00%) 0 (0.00%) 1 (6.67%) 2 (10.0%) 
export2xls(export_table,paste0(home,"/docs/Descriptives.xlsx"))

## compare Binary Spatial Controls 
res = compareGroups(Control_spatial_binary~., data = samplemeta[,variables], max.ylev = 10)
Warning in chisq.test(xx, correct = FALSE): Chi-squared approximation may be incorrect

Warning in chisq.test(xx, correct = FALSE): Chi-squared approximation may be incorrect

Warning in chisq.test(xx, correct = FALSE): Chi-squared approximation may be incorrect

Warning in chisq.test(xx, correct = FALSE): Chi-squared approximation may be incorrect

Warning in chisq.test(xx, correct = FALSE): Chi-squared approximation may be incorrect

Warning in chisq.test(xx, correct = FALSE): Chi-squared approximation may be incorrect
#summary(res)
export_table <- createTable(res)
options(width = 10000)
export2md(export_table)
Summary descriptives table by groups of `Control_spatial_binary’
FALSE TRUE p.overall
N=505 N=44
Celltype: 0.080
Microglia_1 338 (66.9%) 33 (75.0%)
Microglia_2 75 (14.9%) 11 (25.0%)
Microglia_3 57 (11.3%) 0 (0.00%)
Microglia_4 8 (1.58%) 0 (0.00%)
Granulocytes 4 (0.79%) 0 (0.00%)
Microglia_5 15 (2.97%) 0 (0.00%)
Mono/Mac 8 (1.58%) 0 (0.00%)
Sex: 0.046
f 118 (23.4%) 4 (9.09%)
m 387 (76.6%) 40 (90.9%)
Age: 0.019
37 weeks 316 (62.6%) 27 (61.4%)
38 weeks 75 (14.9%) 1 (2.27%)
40 weeks 114 (22.6%) 16 (36.4%)
Genotype_corr: <0.001
WT 274 (54.3%) 9 (20.5%)
APPPS1+ 231 (45.7%) 35 (79.5%)
Treatment: 0.795
Ctrl 375 (74.3%) 34 (77.3%)
Stroke 130 (25.7%) 10 (22.7%)
condition: <0.001
APPPS1+_Ctrl 178 (35.2%) 29 (65.9%)
APPPS1+_Stroke 53 (10.5%) 6 (13.6%)
WT_Ctrl 197 (39.0%) 5 (11.4%)
WT_Stroke 77 (15.2%) 4 (9.09%)
Phase: 0.215
G1 217 (43.0%) 13 (29.5%)
G2M 123 (24.4%) 14 (31.8%)
S 165 (32.7%) 17 (38.6%)
Stroke_spatial_binary: <0.001
FALSE 371 (73.5%) 8 (18.2%)
TRUE 134 (26.5%) 36 (81.8%)
Brain_region: 0.795
Cortex 375 (74.3%) 34 (77.3%)
Lesion 130 (25.7%) 10 (22.7%)
nCount_RNA 161418 (112945) 214465 (106144) 0.003
pseudoaligned_reads 161849 (113014) 214704 (106166) 0.003
percent.mt 1.76 (1.27) 3.14 (2.04) <0.001
percent.rp 2.05 (1.34) 4.20 (2.00) <0.001
Mouse_ID: .
256#1022 58 (11.5%) 2 (4.55%)
256#1023 64 (12.7%) 2 (4.55%)
364#469 75 (14.9%) 1 (2.27%)
386 19 (3.76%) 1 (2.27%)
387 10 (1.98%) 3 (6.82%)
388 85 (16.8%) 12 (27.3%)
409 58 (11.5%) 3 (6.82%)
457 93 (18.4%) 17 (38.6%)
461 43 (8.51%) 3 (6.82%) 
export2xls(export_table,paste0(home,"/docs/Descriptives_sp_Control.xlsx"))

## compare Binary Spatial  Stroke
res = compareGroups(Stroke_spatial_binary~., data = samplemeta[,variables], max.ylev = 10)
#summary(res)
export_table <- createTable(res)
options(width = 10000)
export2md(export_table)
Summary descriptives table by groups of `Stroke_spatial_binary’
FALSE TRUE p.overall
N=379 N=170
Celltype: .
Microglia_1 240 (63.3%) 131 (77.1%)
Microglia_2 55 (14.5%) 31 (18.2%)
Microglia_3 53 (14.0%) 4 (2.35%)
Microglia_4 8 (2.11%) 0 (0.00%)
Granulocytes 4 (1.06%) 0 (0.00%)
Microglia_5 11 (2.90%) 4 (2.35%)
Mono/Mac 8 (2.11%) 0 (0.00%)
Sex: <0.001
f 103 (27.2%) 19 (11.2%)
m 276 (72.8%) 151 (88.8%)
Age: <0.001
37 weeks 242 (63.9%) 101 (59.4%)
38 weeks 76 (20.1%) 0 (0.00%)
40 weeks 61 (16.1%) 69 (40.6%)
Genotype_corr: <0.001
WT 249 (65.7%) 34 (20.0%)
APPPS1+ 130 (34.3%) 136 (80.0%)
Treatment: 0.003
Ctrl 297 (78.4%) 112 (65.9%)
Stroke 82 (21.6%) 58 (34.1%)
condition: <0.001
APPPS1+_Ctrl 97 (25.6%) 110 (64.7%)
APPPS1+_Stroke 33 (8.71%) 26 (15.3%)
WT_Ctrl 200 (52.8%) 2 (1.18%)
WT_Stroke 49 (12.9%) 32 (18.8%)
Phase: 0.573
G1 155 (40.9%) 75 (44.1%)
G2M 93 (24.5%) 44 (25.9%)
S 131 (34.6%) 51 (30.0%)
Control_spatial_binary: <0.001
FALSE 371 (97.9%) 134 (78.8%)
TRUE 8 (2.11%) 36 (21.2%)
Brain_region: 0.003
Cortex 297 (78.4%) 112 (65.9%)
Lesion 82 (21.6%) 58 (34.1%)
nCount_RNA 161384 (122017) 175222 (90282) 0.139
pseudoaligned_reads 161834 (122084) 175563 (90328) 0.143
percent.mt 1.54 (1.03) 2.62 (1.78) <0.001
percent.rp 1.59 (1.06) 3.62 (1.46) <0.001
Mouse_ID: .
256#1022 58 (15.3%) 2 (1.18%)
256#1023 66 (17.4%) 0 (0.00%)
364#469 76 (20.1%) 0 (0.00%)
386 8 (2.11%) 12 (7.06%)
387 6 (1.58%) 7 (4.12%)
388 47 (12.4%) 50 (29.4%)
409 41 (10.8%) 20 (11.8%)
457 50 (13.2%) 60 (35.3%)
461 27 (7.12%) 19 (11.2%) 
export2xls(export_table,paste0(home,"/docs/Descriptives_sp_Stroke.xlsx"))

download data as excel file here

Preprocessing

include only genes that are different across samples (5 974 excluded) include only genes with more than 10 reads in at least 10 cells in at least one Celltype (21 733 transcripts excluded) exclude all cells with less than 10000 reads (none excluded)

649 cells and 11 623 transcripts analyzed

#get normalized counts 
# question to Desiree hat the Seurat object been initialized with normalized data?
counts <- samples.integrated@assays$RNA@counts %>% as.data.frame()

# drop no variance data and sort by samplemeta
counts <- counts[apply(counts,1, sd) > 0, rownames(samplemeta)]

# drop genes with low detection rate (more than 5 counts per cell)
counts_per_celltype=apply(counts, 1, function(x){tapply(x, samplemeta$Celltype, function(z){sum(z>10,na.rm=T)})}) %>% as.data.frame()

# keep RNAs with at least 10 cells with good expression
idxr=which(colSums(counts_per_celltype)>10)

idxc=which(colSums(counts)>10000)

counts = counts[idxr,idxc]
samplemeta=samplemeta[colnames(counts), ]

Overview tables

samplemeta$corrGenotype_Treatment<-
  paste0(samplemeta$Genotype_corr," ", samplemeta$Treatment)

sumstat = counts %>% t() %>%  
  aggregate(by=samplemeta["corrGenotype_Treatment"], mean) %>% 
  t() %>% as.data.frame()

names(sumstat) <- paste0("mean",sumstat[1,])
sumstat<-sumstat[-1,]

sumstatCelltype = counts %>% t() %>%  
  aggregate(by=samplemeta["Celltype"], mean) %>% 
  t() %>% as.data.frame()

names(sumstatCelltype) <- paste0("mean",sumstatCelltype[1,])
sumstatCelltype<-sumstatCelltype[-1,]

sumstat_spStroke = counts %>% t() %>%  
  aggregate(by=samplemeta["Stroke_spatial_binary"], mean) %>% 
  t() %>% as.data.frame()

names(sumstat_spStroke) <- paste0("mean_spStrk",sumstat_spStroke[1,])
sumstat_spStroke<-sumstat_spStroke[-1,]


sumstat_spCtrl = counts %>% t() %>%  
  aggregate(by=samplemeta["Control_spatial_binary"], mean) %>% 
  t() %>% as.data.frame()

names(sumstat_spCtrl) <- paste0("mean_spStrk",sumstat_spCtrl[1,])
sumstat_spCtrl<-sumstat_spCtrl[-1,]



Geno.Treatment_kruskal=apply(counts,1, function(x){
  res = kruskal.test(unlist(x)~samplemeta$corrGenotype_Treatment) %>% unlist()
  return(as.numeric(res["p.value"]))
})

Geno.Treatment_kruskal_fdr = p.adjust(Geno.Treatment_kruskal, method = "fdr")

Celltype_kruskal=apply(counts,1, function(x){
   res = kruskal.test(unlist(x)~samplemeta$Celltype) %>% unlist()
  return(as.numeric(res["p.value"]))
})

Celltype_kruskal_fdr = p.adjust(Celltype_kruskal, method = "fdr")
                                
                                
Control_spatial_binary_kruskal=apply(counts,1, function(x){
  res = kruskal.test(unlist(x)~samplemeta$Control_spatial_binary) %>% unlist()
  return(as.numeric(res["p.value"]))
})
Control_spatial_binary_kruskal_fdr = p.adjust(Control_spatial_binary_kruskal, method = "fdr")

Stroke_spatial_binary_kruskal=apply(counts,1, function(x){
  res = kruskal.test(unlist(x)~samplemeta$Stroke_spatial_binary) %>% unlist()
  return(as.numeric(res["p.value"]))
})

Stroke_spatial_binary_kruskal_fdr = p.adjust(Stroke_spatial_binary_kruskal, method = "fdr")


summarytab=cbind(sumstat, Geno.Treatment_kruskal, Geno.Treatment_kruskal_fdr , 
                 sumstatCelltype, Celltype_kruskal, Celltype_kruskal_fdr, 
                 sumstat_spCtrl, Control_spatial_binary_kruskal, Control_spatial_binary_kruskal_fdr,
                 sumstat_spStroke,Stroke_spatial_binary_kruskal,Stroke_spatial_binary_kruskal_fdr)

summarytab %>% display_tab()
Warning in instance$preRenderHook(instance): It seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.html

Hierarchical clustering

to check where the variance in the data comes from

log2_cpm = log2(counts+1)

varsset=apply(log2_cpm, 1, var)

cpm.sel.trans = t(log2_cpm[order(varsset,decreasing = T)[1:1000],])

distance = dist(cpm.sel.trans)

sampleDistMatrix <- as.matrix(distance)

#colors for plotting heatmap
colors <- rev(colorRampPalette(brewer.pal(9, "Spectral"))(255))
colors=jetcolors(255)
colors=viridis(255)

cellcol = Dark8[1:nlevels(samplemeta$Celltype)]
names(cellcol) = levels(samplemeta$Celltype)

Stroke_sp = Dark8[1:nlevels(samplemeta$Stroke_spatial_binary)]
names(Stroke_sp) = levels(samplemeta$Stroke_spatial_binary)

Ctrl_sp = Dark8[1:nlevels(samplemeta$Control_spatial_binary)]
names(Ctrl_sp) = levels(samplemeta$Control_spatial_binary)

genotypecol = brewer.pal(4,"Accent")[c(1:nlevels(samplemeta$Genotype_corr))]
names(genotypecol) = levels(samplemeta$Genotype_corr)

strokecol = brewer.pal(5,"Set2")[1:nlevels(samplemeta$Treatment)+2]
names(strokecol) = levels(samplemeta$Treatment)

mousecol = brewer.pal(9,"Set1")[1:nlevels(samplemeta$Mouse_ID)]
names(mousecol) = levels(samplemeta$Mouse_ID)

braincol = brewer.pal(3,"Set2")[1:nlevels(samplemeta$Brain_region)]
names(braincol) = levels(samplemeta$Brain_region)

samplemeta$Age = as.factor(samplemeta$Age)
Agecol = colorRampPalette(c("dodgerblue",
                          "dodgerblue4"))(nlevels(samplemeta$Age))[1:nlevels(samplemeta$Age)]
names(Agecol) = levels(samplemeta$Age)

ann_colors = list(
  Genotype_corr = genotypecol, 
  Mouse_ID = mousecol,
  Brain_region = braincol,
  Celltype=cellcol,
  Treatment=strokecol, 
  Age=Agecol, 
  Control_spatial_binary = Ctrl_sp,
  Stroke_spatial_binary =Stroke_sp)

labels = samplemeta[,c("Genotype_corr","Mouse_ID", "Age", "Brain_region", "Celltype", 
                       "Control_spatial_binary","Stroke_spatial_binary","Treatment")] %>%  
  mutate_all(as.character) %>% as.data.frame()
rownames(labels)=rownames(samplemeta)

pheatmap(sampleDistMatrix,
         clustering_distance_rows = distance,
         clustering_distance_cols = distance,
         clustering_method = "ward.D2",
         scale ="none",
         legend = F,
         show_rownames=F, show_colnames = F,
         border_color = NA, 
         annotation_row = labels,
         annotation_col = labels,
         annotation_colors = ann_colors,
         col = colors, 
         main = "D62 top1000 Distances normalized log2 counts")

Version Author Date
d5c8125 Andreas Chiocchetti 2023-07-14
982b9f6 Andreas Chiocchetti 2023-02-20
eaa373d Andreas Chiocchetti 2023-02-16
cfbfcf6 achiocch 2022-06-22
51b35c8 achiocch 2022-04-26
8e44cfc achiocch 2022-04-22
55e7738 achiocch 2022-03-30
f21fdf7 achiocch 2022-03-30
858d646 achiocch 2022-02-03

Version Author Date
d5c8125 Andreas Chiocchetti 2023-07-14
eaa373d Andreas Chiocchetti 2023-02-16
cfbfcf6 achiocch 2022-06-22

Genotype X Treatment Statistical modelling

the logFC should not be interpreted on its own without the specific post hoc tests (see app below to check the individual genes) In any case the value here would correspond to the APPPS1+ with Stroke against all others.

Heatmaps show standardized deviations from mean across all cells, trimmed to 2 standard deviations ( i.e. values above 2 SD are set to 2 SD). this allows to see more subtle changes better.

getres=function(Celltype="Specify", 
                Hypothesis="~1+Sex+Genotype_corr*Treatment", 
                Target="Genotype_corrAPPPS1+:TreatmentStroke",
                Randomeffect=NULL){
  
  res= comparison_rand(designform =Hypothesis, 
                       randomeffect = Randomeffect, 
                       Samples = samplemeta$Celltype==Celltype, 
                       log_cpm = log2_cpm, 
                       samplesdata = samplemeta,
                       target=Target)
  
  
  
  labels = samplemeta[samplemeta$Celltype==Celltype,c("Treatment","Age", "Genotype_corr","Mouse_ID", "Brain_region", "Celltype","Control_spatial_binary", "Stroke_spatial_binary")] %>%  
    mutate_all(as.character) %>% as.data.frame()

  
  labels=labels %>% arrange(Treatment,Genotype_corr)
  plotdata=log2_cpm[res$adj.P.Val<0.05,rownames(labels)]
  plotdata = apply(plotdata, 2, trimmed_scaled)
  
  
  
  pheatmap(plotdata,
           #clustering_method = "ward.D2",
           cluster_cols = F,
           cluster_rows  = T,
           scale ="row",
           show_rownames=F, show_colnames = F,
           legend=T,
           border_color = NA, 
           #annotation_row = labels,
           annotation_col = labels,
           annotation_colors = ann_colors,
           col = colors, 
           breaks = seq(-2,2, length.out=254),
           main = "D62 Distances normalized log2 counts")
  return(res)
  
}

getresspatio=function(Celltype="Specify", 
                Hypothesis="~1+Sex+Genotype_corr*Treatment", 
                Target="Genotype_corrAPPPS1+:TreatmentStroke",
                Randomeffect=NULL){
  
  idx=NA
  if(Celltype =="Stroke_spatial_binary"){
    idx=samplemeta$Stroke_spatial_binary==T
  
  }
  if(Celltype =="Control_spatial_binary"){
    idx=samplemeta$Control_spatial_binary==T
    }
  
  idx[is.na(idx)]=F

  res=comparison_rand(designform =Hypothesis, 
                       randomeffect = Randomeffect, 
                       Samples = idx, 
                       log_cpm = log2_cpm, 
                       samplesdata = samplemeta,
                       target=Target)
  
  
  
  labels = samplemeta[idx,c("Treatment","Age", "Genotype_corr","Mouse_ID", "Brain_region", "Celltype", "Control_spatial_binary", "Stroke_spatial_binary")] %>%  
    mutate_all(as.character) %>% as.data.frame()

  
  labels=labels %>% arrange(Treatment,Genotype_corr)
  
  plotdata=log2_cpm[res$adj.P.Val<0.05,rownames(labels)]
  plotdata = apply(plotdata, 2, trimmed_scaled)
  
  
  
  pheatmap(plotdata,
           #clustering_method = "ward.D2",
           cluster_cols = F,
           cluster_rows  = T,
           scale ="row",
           show_rownames=F, show_colnames = F,
           legend=T,
           border_color = NA, 
           #annotation_row = labels,
           annotation_col = labels,
           annotation_colors = ann_colors,
           col = colors, 
           breaks = seq(-2,2, length.out=254),
           main = "D62 Distances normalized log2 counts")
  
  return(res)
  
}

generate_output=function(CT="Celltype",mode="Celltype", ...){
  respath=paste0(home, "/docs/LMER_",CT,".xlsx")
  
  if(mode == "Spatial"){
    analysis = getresspatio(Celltype = CT)
  }
  if (mode == "Celltype"){
    analysis = getres(Celltype = CT)}
  
  analysis_sig = analysis[analysis$adj.P.Val<0.05,]
  analysis_sig %>% display_tab()
  
  write.xlsx2(analysis_sig, file =respath , sheetName = "significant genes")
  
  ResGO = getGOresults(rownames(analysis_sig),
                       rownames(analysis),
                       "mmusculus")
  
  if(length(ResGO)>0){
    p=gostplot(ResGO)
  } else{
    ResGO=data.frame(result="no significant enrichment identified")
    p="no significant enrichment identified"
  }
  
  write.xlsx2(ResGO$result, file = respath, sheetName = "GO_enrichment", append=T)
  return(list(results=analysis,results_sig=analysis_sig, plot=p))
}

Microglia_1

res_output=generate_output(CT = "Microglia_1")

Version Author Date
d5c8125 Andreas Chiocchetti 2023-07-14
982b9f6 Andreas Chiocchetti 2023-02-20
eaa373d Andreas Chiocchetti 2023-02-16
cfbfcf6 achiocch 2022-06-22
8e44cfc achiocch 2022-04-22
55e7738 achiocch 2022-03-30
87439a3 achiocch 2022-02-21
79dd0b1 achiocch 2022-02-21
858d646 achiocch 2022-02-03 
[1] "no significant GO terms identified"
res_output[["results_sig"]] %>% display_tab()
res_output[["plot"]]
[1] "no significant enrichment identified"

Data sources and their abbreviations are: Gene Ontology (GO or by branch GO:MF, GO:BP, GO:CC)KEGG (KEGG) TRANSFAC (TF) miRTarBase (MIRNA) CORUM (CORUM) Human phenotype ontology (HP) Human Protein Atlas (HPA)

get results table here

Microglia_2

res_output=generate_output("Microglia_2")

Version Author Date
d5c8125 Andreas Chiocchetti 2023-07-14
982b9f6 Andreas Chiocchetti 2023-02-20
eaa373d Andreas Chiocchetti 2023-02-16
cfbfcf6 achiocch 2022-06-22
8e44cfc achiocch 2022-04-22
55e7738 achiocch 2022-03-30
87439a3 achiocch 2022-02-21
79dd0b1 achiocch 2022-02-21
858d646 achiocch 2022-02-03 
[1] "no significant GO terms identified"
res_output[["results_sig"]] %>% display_tab()
res_output[["plot"]]
[1] "no significant enrichment identified"

Data sources and their abbreviations are: Gene Ontology (GO or by branch GO:MF, GO:BP, GO:CC)KEGG (KEGG) TRANSFAC (TF) miRTarBase (MIRNA) CORUM (CORUM) Human phenotype ontology (HP) Human Protein Atlas (HPA)

get results table here

Microglia_3

idx=samplemeta$Celltype=="Microglia_3"
table(samplemeta$Treatment[idx], samplemeta$Genotype_corr[idx])
        
         WT APPPS1+
  Ctrl   42      10
  Stroke  0       5
table(samplemeta$Treatment[idx], samplemeta$Genotype_corr[idx]) %>% fisher.test()

    Fisher's Exact Test for Count Data

data:  .
p-value = 0.0007172
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 3.138602      Inf
sample estimates:
odds ratio 
       Inf

Data sources and their abbreviations are: Gene Ontology (GO or by branch GO:MF, GO:BP, GO:CC)KEGG (KEGG) TRANSFAC (TF) miRTarBase (MIRNA) CORUM (CORUM) Human phenotype ontology (HP) Human Protein Atlas (HPA)

get results table here

Microglia_4

calculation not converging as Celltype not identified in all conditions however also rare in other conditions, no significant difference

idx=samplemeta$Celltype=="Microglia_4"
table(samplemeta$Treatment[idx], samplemeta$Genotype_corr[idx])
        
         WT APPPS1+
  Ctrl    7       0
  Stroke  0       0
table(samplemeta$Treatment[idx], samplemeta$Genotype_corr[idx]) %>% fisher.test()

    Fisher's Exact Test for Count Data

data:  .
p-value = 1
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
   0 Inf
sample estimates:
odds ratio 
         0 
#res_output=generate_output("Microglia_4")
#res_output[["results_sig"]] %>% display_tab()
#res_output[["plot"]]

get results table here

Microglia_5

res_output=generate_output("Microglia_5")

Version Author Date
d5c8125 Andreas Chiocchetti 2023-07-14
982b9f6 Andreas Chiocchetti 2023-02-20
eaa373d Andreas Chiocchetti 2023-02-16
cfbfcf6 achiocch 2022-06-22
8e44cfc achiocch 2022-04-22
55e7738 achiocch 2022-03-30
87439a3 achiocch 2022-02-21
79dd0b1 achiocch 2022-02-21
858d646 achiocch 2022-02-03 
[1] "no significant GO terms identified"
res_output[["results_sig"]] %>% display_tab()
res_output[["plot"]]
[1] "no significant enrichment identified"

Data sources and their abbreviations are: Gene Ontology (GO or by branch GO:MF, GO:BP, GO:CC)KEGG (KEGG) TRANSFAC (TF) miRTarBase (MIRNA) CORUM (CORUM) Human phenotype ontology (HP) Human Protein Atlas (HPA) get results table here

Granulozytes

calculation not converging as Celltype not identified in all conditions however also rare in other conditions, no significant difference

idx=samplemeta$Celltype=="Granulocytes"
table(samplemeta$Treatment[idx], samplemeta$Genotype_corr[idx])
        
         WT APPPS1+
  Ctrl    4       0
  Stroke  0       0
table(samplemeta$Treatment[idx], samplemeta$Genotype_corr[idx]) %>% fisher.test()

    Fisher's Exact Test for Count Data

data:  .
p-value = 1
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
   0 Inf
sample estimates:
odds ratio 
         0 
# res_output=generate_output("Granulocytes ")
# res_output[["results_sig"]] %>% display_tab()
# res_output[["plot"]]

get results table here

Mono/Mac

idx=samplemeta$Celltype=="Mono/Mac"
table(samplemeta$Treatment[idx], samplemeta$Genotype_corr[idx])
        
         WT APPPS1+
  Ctrl    7       1
  Stroke  8       4
table(samplemeta$Treatment[idx], samplemeta$Genotype_corr[idx]) %>% fisher.test()

    Fisher's Exact Test for Count Data

data:  .
p-value = 0.6027
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
   0.2433699 196.3003974
sample estimates:
odds ratio 
  3.302108 
# res_output=generate_output("T/NK")
# res_output[["results_sig"]] %>% display_tab()
# res_output[["plot"]]

get results table here

Control Spacial

res_output=generate_output(CT = "Control_spatial_binary", mode="Spatial")

Version Author Date
d5c8125 Andreas Chiocchetti 2023-07-14 
[1] "no significant GO terms identified"
res_output[["results_sig"]] %>% display_tab()
res_output[["plot"]]
[1] "no significant enrichment identified"

Strole Spatial

res_output=generate_output("Stroke_spatial_binary", mode="Spatial")

Version Author Date
d5c8125 Andreas Chiocchetti 2023-07-14 
[1] "no significant GO terms identified"
res_output[["results_sig"]] %>% display_tab()
res_output[["plot"]]
[1] "no significant enrichment identified"

Pseudotime analysis microglia only

identify starting cluster

Riddgplot suggests Microglia_4 as the homeostatic one

samples.microglia = samples.integrated[,grepl("Microgl", samples.integrated$Celltype)]

RidgePlot(samples.microglia, features = c("Tmem119", "Cx3cr1", "P2ry12", "P2ry13", "Axl", "Itgax", "Cst7", "Clec7a", "Spp1"), ncol = 2)
Picking joint bandwidth of 0.581
Picking joint bandwidth of 0.247
Picking joint bandwidth of 0.441
Picking joint bandwidth of 0.598
Picking joint bandwidth of 0.431
Picking joint bandwidth of 0.0287
Picking joint bandwidth of 0.601
Picking joint bandwidth of 0.524
Picking joint bandwidth of 5.04e-06

Version Author Date
9d29eff Andreas Chiocchetti 2023-08-11 
samplemeta.microglia = samples.microglia@meta.data
sce.microglia=as.SingleCellExperiment(samples.microglia)
sce.microglia <- slingshot(sce.microglia, clusterLabels = 'Celltype', reducedDim = 'UMAP', start.clus="Microglia_4")

samplemeta.microglia$Pseudotime_1 <- sce.microglia$slingPseudotime_1
samplemeta.microglia$Pseudotime_2 <- sce.microglia$slingPseudotime_2


dataset = data.frame(reducedDim(sce.microglia, "UMAP"), 
                     cluster = as.character(sce.microglia$Celltype), 
                     Stroke_spatial = as.character(sce.microglia$Stroke_spatial_binary))



curves <- slingCurves(sce.microglia, as.df=T)

p = ggplot(dataset, aes(x=UMAP_1, y=UMAP_2))+
  geom_point(aes(col = cluster), size=2)+
  geom_point(data=subset(dataset, dataset$Stroke_spatial==T), aes(fill=Stroke_spatial), pch=21, size=0.7)+
  scale_fill_manual(values="black")
curves = curves %>% arrange(Lineage, Order)

p = p + geom_path(data = curves, aes(group = Lineage), lwd=1.5, lineend = "butt",)+theme_classic()

curves_arrows = curves %>% group_by(Lineage) %>%  mutate(xend = lead(UMAP_1), yend = lead(UMAP_2))

curves_arrows = curves_arrows[curves_arrows$Order %in% seq(10,150, by=floor((150-10)/10)),]

p = p + geom_segment(data = curves_arrows,aes(x = UMAP_1, y = UMAP_2, xend = xend, yend = yend), arrow = arrow(length=unit(0.3, "cm"), type="closed"))
p
Warning: Removed 2 rows containing missing values (`geom_segment()`).

Version Author Date
9d29eff Andreas Chiocchetti 2023-08-11
87fbacb Andreas Chiocchetti 2023-08-07
982b9f6 Andreas Chiocchetti 2023-02-20
eaa373d Andreas Chiocchetti 2023-02-16
5bba014 achiocch 2022-06-27
cfbfcf6 achiocch 2022-06-22 
ggsave(filename = "./docs/Pseudotime_linages.svg", p)
Saving 6 x 6 in image
Warning: Removed 2 rows containing missing values (`geom_segment()`).
ggsave(filename = "./docs/Pseudotime_linages.pdf", p)
Saving 6 x 6 in image
Warning: Removed 2 rows containing missing values (`geom_segment()`).
dataset$pseudotime_1 <- sce.microglia$slingPseudotime_1
dataset$pseudotime_2 <- sce.microglia$slingPseudotime_2
dataset$Genotype_Treatment <- sce.microglia$Genotype_Treatment

metadata = samples.integrated@meta.data %>%  as.data.frame()
metadata$Pseudotime_1=NA
metadata$Pseudotime_2=NA

metadata[rownames(dataset),]$Pseudotime_1 = dataset$pseudotime_1
metadata[rownames(dataset),]$Pseudotime_2 = dataset$pseudotime_2

samples.integrated@meta.data <- metadata
linages = sce.microglia %>% slingLineages()


res = compareGroups(Celltype~., data = metadata[,c(variables, "Pseudotime_1", "Pseudotime_2")], max.ylev = 10)
export_table <- createTable(res)
options(width = 10000)
export2md(export_table)
Summary descriptives table by groups of `Celltype’
Microglia_1 Microglia_2 Microglia_3 Microglia_4 Granulocytes Microglia_5 Mono/Mac p.overall
N=372 N=86 N=57 N=8 N=4 N=15 N=20
Sex: 0.276
f 73 (19.6%) 20 (23.3%) 17 (29.8%) 2 (25.0%) 2 (50.0%) 5 (33.3%) 5 (25.0%)
m 299 (80.4%) 66 (76.7%) 40 (70.2%) 6 (75.0%) 2 (50.0%) 10 (66.7%) 15 (75.0%)
Age: .
37 weeks 226 (60.8%) 56 (65.1%) 41 (71.9%) 6 (75.0%) 2 (50.0%) 7 (46.7%) 11 (55.0%)
38 weeks 44 (11.8%) 8 (9.30%) 13 (22.8%) 2 (25.0%) 2 (50.0%) 4 (26.7%) 3 (15.0%)
40 weeks 102 (27.4%) 22 (25.6%) 3 (5.26%) 0 (0.00%) 0 (0.00%) 4 (26.7%) 6 (30.0%)
Genotype_corr: .
WT 170 (45.7%) 43 (50.0%) 42 (73.7%) 8 (100%) 4 (100%) 8 (53.3%) 15 (75.0%)
APPPS1+ 202 (54.3%) 43 (50.0%) 15 (26.3%) 0 (0.00%) 0 (0.00%) 7 (46.7%) 5 (25.0%)
Treatment: .
Ctrl 272 (73.1%) 56 (65.1%) 52 (91.2%) 8 (100%) 4 (100%) 11 (73.3%) 8 (40.0%)
Stroke 100 (26.9%) 30 (34.9%) 5 (8.77%) 0 (0.00%) 0 (0.00%) 4 (26.7%) 12 (60.0%)
condition: .
APPPS1+_Ctrl 164 (44.1%) 28 (32.6%) 10 (17.5%) 0 (0.00%) 0 (0.00%) 6 (40.0%) 1 (5.00%)
APPPS1+_Stroke 38 (10.2%) 15 (17.4%) 5 (8.77%) 0 (0.00%) 0 (0.00%) 1 (6.67%) 4 (20.0%)
WT_Ctrl 108 (29.0%) 28 (32.6%) 42 (73.7%) 8 (100%) 4 (100%) 5 (33.3%) 7 (35.0%)
WT_Stroke 62 (16.7%) 15 (17.4%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 3 (20.0%) 8 (40.0%)
Phase: .
G1 159 (42.7%) 36 (41.9%) 19 (33.3%) 0 (0.00%) 2 (50.0%) 10 (66.7%) 9 (45.0%)
G2M 95 (25.5%) 18 (20.9%) 17 (29.8%) 2 (25.0%) 2 (50.0%) 2 (13.3%) 6 (30.0%)
S 118 (31.7%) 32 (37.2%) 21 (36.8%) 6 (75.0%) 0 (0.00%) 3 (20.0%) 5 (25.0%)
Control_spatial_binary: 0.080
FALSE 338 (91.1%) 75 (87.2%) 57 (100%) 8 (100%) 4 (100%) 15 (100%) 8 (100%)
TRUE 33 (8.89%) 11 (12.8%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 0 (0.00%)
Stroke_spatial_binary: .
FALSE 240 (64.7%) 55 (64.0%) 53 (93.0%) 8 (100%) 4 (100%) 11 (73.3%) 8 (100%)
TRUE 131 (35.3%) 31 (36.0%) 4 (7.02%) 0 (0.00%) 0 (0.00%) 4 (26.7%) 0 (0.00%)
Brain_region: .
Cortex 272 (73.1%) 56 (65.1%) 52 (91.2%) 8 (100%) 4 (100%) 11 (73.3%) 8 (40.0%)
Lesion 100 (26.9%) 30 (34.9%) 5 (8.77%) 0 (0.00%) 0 (0.00%) 4 (26.7%) 12 (60.0%)
nCount_RNA 158416 (103720) 194411 (106729) 164921 (154435) 217767 (205467) 204215 (86217) 165673 (126013) 162375 (90426) 0.154
pseudoaligned_reads 158947 (104087) 194836 (106825) 165764 (154492) 218138 (205598) 204425 (86225) 166053 (126105) 163005 (90604) 0.160
percent.mt 1.90 (1.35) 2.14 (1.78) 1.44 (1.07) 1.66 (1.04) 0.98 (1.13) 2.01 (1.01) 2.04 (1.27) 0.089
percent.rp 2.21 (1.29) 3.05 (2.17) 1.48 (1.34) 0.98 (0.41) 1.02 (0.50) 2.03 (1.25) 3.16 (2.30) <0.001
Mouse_ID: .
256#1022 33 (8.87%) 9 (10.5%) 13 (22.8%) 2 (25.0%) 0 (0.00%) 1 (6.67%) 2 (10.0%)
256#1023 31 (8.33%) 11 (12.8%) 16 (28.1%) 4 (50.0%) 2 (50.0%) 0 (0.00%) 2 (10.0%)
364#469 44 (11.8%) 8 (9.30%) 13 (22.8%) 2 (25.0%) 2 (50.0%) 4 (26.7%) 3 (15.0%)
386 16 (4.30%) 4 (4.65%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 4 (20.0%)
387 9 (2.42%) 3 (3.49%) 1 (1.75%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 2 (10.0%)
388 77 (20.7%) 15 (17.4%) 2 (3.51%) 0 (0.00%) 0 (0.00%) 4 (26.7%) 0 (0.00%)
409 46 (12.4%) 11 (12.8%) 0 (0.00%) 0 (0.00%) 0 (0.00%) 3 (20.0%) 4 (20.0%)
457 87 (23.4%) 13 (15.1%) 8 (14.0%) 0 (0.00%) 0 (0.00%) 2 (13.3%) 1 (5.00%)
461 29 (7.80%) 12 (14.0%) 4 (7.02%) 0 (0.00%) 0 (0.00%) 1 (6.67%) 2 (10.0%)
Pseudotime_1 10.7 (1.86) 5.17 (1.34) 3.89 (0.17) 0.48 (0.48) . (.) 14.2 (0.29) . (.) <0.001
Pseudotime_2 8.08 (1.90) 5.46 (1.18) 8.94 (1.36) 0.48 (0.48) . (.) . (.) . (.) <0.001 
export2xls(export_table,paste0(home,"/docs/Descriptives_withPseudotime.xlsx"))


p1 <-ggplot(dataset, aes(x=pseudotime_1))+
  geom_density(aes(group=Genotype_Treatment, col=Genotype_Treatment), lwd=2)+
  theme_classic()

p = ggplot(dataset, aes(x=UMAP_1, y=UMAP_2))+
  geom_point(aes(col = pseudotime_1), size=4)+
  geom_point(data=subset(dataset, dataset$Stroke_spatial==T), aes(fill=Stroke_spatial), pch=21)+
  scale_fill_manual(values="black")+scale_color_viridis()+theme_classic()
p = p + geom_path(data = curves[curves$Lineage==1,] %>% arrange(Order), aes(group = Lineage), lwd=2)+theme_classic()
p = p + geom_segment(data = curves_arrows[curves_arrows$Lineage==1,],aes(x = UMAP_1, y = UMAP_2, xend = xend, yend = yend), arrow = arrow(length=unit(0.3, "cm"), type="closed"))

plt = ggarrange(p1,p)
Warning: Removed 67 rows containing non-finite values (`stat_density()`).
Warning: Removed 1 rows containing missing values (`geom_segment()`).
plt

Version Author Date
9d29eff Andreas Chiocchetti 2023-08-11
87fbacb Andreas Chiocchetti 2023-08-07 
ggsave(filename = "./docs/Pseudotime_density_Lin1.svg", plt)
Saving 12 x 6 in image
ggsave(filename = "./docs/Pseudotime_denisty_Lin1.pdf", plt)
Saving 12 x 6 in image
p1 <-ggplot(dataset, aes(x=pseudotime_2))+
  geom_density(aes(group=Genotype_Treatment, col=Genotype_Treatment), lwd=2)+
  theme_classic()

p = ggplot(dataset, aes(x=UMAP_1, y=UMAP_2))+
  geom_point(aes(col = pseudotime_2), size=4)+
  geom_point(data=subset(dataset, dataset$Stroke_spatial==T), aes(fill=Stroke_spatial), pch=21)+
  scale_fill_manual(values="black")+scale_color_viridis()+theme_classic()
# manually selected linage 3 here as this fits best to peudotimes values
p = p + geom_path(data = subset(curves, Lineage==2) %>% arrange(Order), aes(group = Lineage), lwd=2)+theme_classic()
p = p + geom_segment(data = curves_arrows[curves_arrows$Lineage==2,],aes(x = UMAP_1, y = UMAP_2, xend = xend, yend = yend), arrow = arrow(length=unit(0.3, "cm"), type="closed"))

plt = ggarrange(p1,p)
Warning: Removed 381 rows containing non-finite values (`stat_density()`).
Warning: Removed 1 rows containing missing values (`geom_segment()`).
plt

Version Author Date
9d29eff Andreas Chiocchetti 2023-08-11
87fbacb Andreas Chiocchetti 2023-08-07 
ggsave(filename = "./docs/Pseudotime_density_Lin2.svg", plt)
Saving 8 x 4 in image
ggsave(filename = "./docs/Pseudotime_denisty_Lin2.pdf", plt)
Saving 8 x 4 in image

Pseudotime Lineage 1 association

GOplot <- function(gostresult)
  if(is.null(gostresult)){
    ggplot() +                      # Draw ggplot2 plot with text only
      annotate("text",
               x = 1,
               y = 1,
               size = 8,
               label = "no significant GO-terms associated") + 
      theme_void()
  } else {
    gostplot(gostresult)}


sce_Pseudotime1 = samples.integrated[,samples.integrated$Celltype %in% linages$Lineage1]

res = comparison_rand(designform ="~1+Sex+Pseudotime_1", 
                Samples = c(), 
                log_cpm = log1p(as.data.frame(sce_Pseudotime1@assays$RNA@counts)), 
                samplesdata = as.data.frame(sce_Pseudotime1@meta.data), 
                randomeffect = c(), 
                target = "Pseudotime_1")
[1] "Sample length is 0 all samples included"
Warning: Zero sample variances detected, have been offset away from zero
analysis_sig <- res[res$adj.P.Val<0.05,]
analysis_sig %>% display_tab()
write.xlsx2(res, file ="./docs/DEX_Pseudotime_Lineage1.xlsx" , sheetName = "DEX genes Lin1")



  
ResGO = getGOresults(rownames(analysis_sig),
                       rownames(res),
                       "mmusculus")

GOplot(ResGO)
ResGO$result %>% display_tab()
sce_Pseudotime1@active.ident <- factor(x = sce_Pseudotime1@active.ident, levels = linages$Lineage1)

plt = DoHeatmap(sce_Pseudotime1, features = rownames(analysis_sig[order(analysis_sig$logFC),]))
Warning in DoHeatmap(sce_Pseudotime1, features = rownames(analysis_sig[order(analysis_sig$logFC), : The following features were omitted as they were not found in the scale.data slot for the SCT assay: Gm9644, Gm17661, Gm11760, Adh1, Gm15947, Tprg, Gm15168, Olfr732, Olig3, Gm10745, Gm44901, Gm8454, Gabrg1
plt

Version Author Date
9d29eff Andreas Chiocchetti 2023-08-11 
ggsave(filename = "./docs/Pseudotime_Lin1_Marker.svg", plt)
Saving 8 x 16 in image
ggsave(filename = "./docs/Pseudotime_Lin1_Marker.pdf", plt)
Saving 8 x 16 in image

Pseudotime Lineage 2 association

sce_Pseudotime2 = samples.integrated[,samples.integrated$Celltype %in% linages$Lineage2]

res = comparison_rand(designform ="~1+Sex+Pseudotime_2", 
                Samples = c(), 
                log_cpm = log1p(as.data.frame(sce_Pseudotime2@assays$RNA@counts)), 
                samplesdata = as.data.frame(sce_Pseudotime2@meta.data), 
                randomeffect = c(), 
                target = "Pseudotime_2")
[1] "Sample length is 0 all samples included"
Warning: Zero sample variances detected, have been offset away from zero
analysis_sig <- res[res$adj.P.Val<0.05,]
analysis_sig %>% display_tab()
write.xlsx2(res, file ="./docs/DEX_Pseudotime_Lineage2.xlsx" , sheetName = "DEX genes Lin1")
  
ResGO = getGOresults(rownames(analysis_sig),
                       rownames(res),
                       "mmusculus")

GOplot(ResGO)
ResGO$result %>% display_tab()
sce_Pseudotime2@active.ident <- factor(x = sce_Pseudotime2@active.ident, levels = linages$Lineage2)

plt = DoHeatmap(sce_Pseudotime2, features = rownames(analysis_sig[order(analysis_sig$logFC),]))
Warning in DoHeatmap(sce_Pseudotime2, features = rownames(analysis_sig[order(analysis_sig$logFC), : The following features were omitted as they were not found in the scale.data slot for the SCT assay: AC117656.2, AC135637.1, 9530053A07Rik, Vmn1r78, Olfr633, Ncapg, Olfr689, Gm20438, Tmem207, Gm12763, Gm37752, Sorcs1, 2700029L08Rik, Gm37970, Slc22a8, Gm7337, Cntnap5b, Gm29250, Olfr1406, Gm37087, Vmn1r217, Slc6a3, Gm26052, AC166056.1, Olfr721-ps1, Gm25100, Gm24975, CAAA01205117.2, Olfr307, Olfr1240, Gm37674, Gm28019, Gm26463, Gm44520, Gm8929, AC120394.2, Itgal, Cdh7, Mug-ps1, Pira2, Runx1t1, 4921534H16Rik, Spink5, Gm34004, Slc2a2, Glp2r, Pcdhb5, Il12b, Speer4cos, 5031434C07Rik, Gm43166, Gm30332, Olfr1, 1700121L16Rik, Efcab6, Gm22008, Gm22587, Olfr288, 4930445B16Rik, D630024D03Rik, Mro, Gm1979, Gm12315, Ciita, Olfr1500, Gm37368, Gm7950, Gm25789, Gm14291, Gm12657, 8430408G22Rik, Gm17661, Gm11760, Gm38093, Gm6634, Gm22852, Gm36388, Gm43690, Gm29170, Tmc8
plt

Version Author Date
9d29eff Andreas Chiocchetti 2023-08-11 
ggsave(filename = "./docs/Pseudotime_Lin2_Marker.svg", plt)
Saving 8 x 16 in image
ggsave(filename = "./docs/Pseudotime_Lin2_Marker.pdf", plt)
Saving 8 x 16 in image

Gene expression plot


sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.1 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                  LC_TIME=en_US.UTF-8           LC_COLLATE=en_US.UTF-8        LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8       LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8           LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8      LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] gprofiler2_0.2.1            lm.beta_1.7-1               pheatmap_1.0.12             RColorBrewer_1.1-3          kableExtra_1.3.4            ggpubr_0.6.0                slingshot_2.6.0             TrajectoryUtils_1.6.0       SingleCellExperiment_1.20.0 princurve_2.1.6             DT_0.27                     viridis_0.6.2               viridisLite_0.4.1           xlsx_0.6.5                  brms_2.18.0                 Rcpp_1.0.10                 compareGroups_4.6.0         data.table_1.14.6           SingleR_2.0.0               SeuratObject_4.1.3          Seurat_4.3.0                forcats_1.0.0               stringr_1.5.0               dplyr_1.1.0                 purrr_1.0.1                 readr_2.1.4                 tidyr_1.3.0                 tibble_3.1.8                ggplot2_3.4.1               tidyverse_1.3.2             DESeq2_1.38.3               SummarizedExperiment_1.28.0 Biobase_2.58.0              MatrixGenerics_1.10.0       matrixStats_0.63.0          GenomicRanges_1.50.2        GenomeInfoDb_1.34.9         IRanges_2.32.0              S4Vectors_0.36.1            BiocGenerics_0.44.0         limma_3.54.1                workflowr_1.7.0            

loaded via a namespace (and not attached):
  [1] rsvd_1.0.5                ica_1.0-3                 svglite_2.1.1             ps_1.7.2                  lmtest_0.9-40             rprojroot_2.0.3           crayon_1.5.2              MASS_7.3-58.2             nlme_3.1-162              backports_1.4.1           posterior_1.3.1           reprex_2.0.2              colourpicker_1.2.0        rlang_1.0.6               XVector_0.38.0            ROCR_1.0-11               readxl_1.4.2              irlba_2.3.5.1             callr_3.7.3               flextable_0.8.5           BiocParallel_1.32.5       bit64_4.0.5               glue_1.6.2                loo_2.5.1                 sctransform_0.3.5         vipor_0.4.5               rstan_2.21.8              parallel_4.2.2            processx_3.8.0            spatstat.sparse_3.0-0     AnnotationDbi_1.60.0      spatstat.geom_3.0-6       haven_2.5.1               tidyselect_1.2.0          fitdistrplus_1.1-8        XML_3.99-0.13             zoo_1.8-11                packrat_0.9.0             distributional_0.3.1      chron_2.3-59              xtable_1.8-4              magrittr_2.0.3            evaluate_0.20             gdtools_0.3.0             cli_3.6.0                 zlibbioc_1.44.0           rstudioapi_0.14           miniUI_0.1.1.1            sp_1.6-0                  whisker_0.4.1             bslib_0.4.2               shinystan_2.6.0           shiny_1.7.4               BiocSingular_1.14.0       xfun_0.37                 askpass_1.1               inline_0.3.19             pkgbuild_1.4.0            cluster_2.1.4             bridgesampling_1.1-2      gfonts_0.2.0              KEGGREST_1.38.0           Brobdingnag_1.2-9         ggrepel_0.9.3             threejs_0.3.3             listenv_0.9.0             xlsxjars_0.6.1            Biostrings_2.66.0         png_0.1-8                 future_1.31.0             withr_2.5.0               bitops_1.0-7              plyr_1.8.8                cellranger_1.1.0          coda_0.19-4               pillar_1.8.1              RcppParallel_5.1.6        cachem_1.0.6              fs_1.6.1                  DelayedMatrixStats_1.20.0 xts_0.12.2                vctrs_0.5.2               ellipsis_0.3.2            generics_0.1.3            dygraphs_1.1.1.6          tools_4.2.2               beeswarm_0.4.0            munsell_0.5.0             DelayedArray_0.24.0       fastmap_1.1.0             compiler_4.2.2            abind_1.4-5               httpuv_1.6.9              plotly_4.10.1             rJava_1.0-6               GenomeInfoDbData_1.2.9    gridExtra_2.3             lattice_0.20-45           deldir_1.0-6              utf8_1.2.3                later_1.3.0               jsonlite_1.8.4            scales_1.2.1              ScaledMatrix_1.6.0        carData_3.0-5             pbapply_1.7-0             sparseMatrixStats_1.10.0  lazyeval_0.2.2            promises_1.2.0.1          car_3.1-1                 goftest_1.2-3             spatstat.utils_3.0-1      reticulate_1.28           checkmate_2.1.0           rmarkdown_2.20            cowplot_1.1.1             textshaping_0.3.6         webshot_0.5.4             Rtsne_0.16                uwot_0.1.14               igraph_1.4.0              rsconnect_0.8.29          survival_3.5-3            yaml_2.3.7                systemfonts_1.0.4         bayesplot_1.10.0          htmltools_0.5.4           rstantools_2.2.0          memoise_2.0.1             locfit_1.5-9.7            digest_0.6.31             assertthat_0.2.1          mime_0.12                 RSQLite_2.2.20            future.apply_1.10.0       blob_1.2.3                ragg_1.2.5                labeling_0.4.2            shinythemes_1.2.0         splines_4.2.2             googledrive_2.0.0         RCurl_1.98-1.10           broom_1.0.3               hms_1.1.2                 modelr_0.1.10             colorspace_2.1-0          base64enc_0.1-3           BiocManager_1.30.19       ggbeeswarm_0.7.1          ggrastr_1.0.1             nnet_7.3-18               sass_0.4.5                RANN_2.6.1                mvtnorm_1.1-3             fansi_1.0.4               tzdb_0.3.0                truncnorm_1.0-8           parallelly_1.34.0         R6_2.5.1                  grid_4.2.2                crul_1.3                  ggridges_0.5.4            lifecycle_1.0.3           StanHeaders_2.21.0-7      zip_2.2.2                 writexl_1.4.2             curl_5.0.0                ggsignif_0.6.4            googlesheets4_1.0.1       leiden_0.4.3              jquerylib_0.1.4           Matrix_1.5-3              RcppAnnoy_0.0.20          spatstat.explore_3.0-6    htmlwidgets_1.6.1         officer_0.5.2             beachmat_2.14.0           polyclip_1.10-4           markdown_1.5              crosstalk_1.2.0           timechange_0.2.0          rvest_1.0.3               globals_0.16.2            openssl_2.0.5             patchwork_1.1.2           spatstat.random_3.1-3     tensorA_0.36.2            progressr_0.13.0          codetools_0.2-19          lubridate_1.9.2           gtools_3.9.4              getPass_0.2-2             prettyunits_1.1.1         dbplyr_2.3.0              gtable_0.3.1              DBI_1.1.3                 git2r_0.31.0              highr_0.10                tensor_1.5                httr_1.4.4                KernSmooth_2.23-20        stringi_1.7.12            reshape2_1.4.4            farver_2.1.1              uuid_1.1-0                annotate_1.76.0           mice_3.15.0               xml2_1.3.3                shinyjs_2.1.0             geneplotter_1.76.0        scattermore_0.8           bit_4.0.5                 spatstat.data_3.0-0       pkgconfig_2.0.3           gargle_1.3.0              rstatix_0.7.2             HardyWeinberg_1.7.5       Rsolnp_1.16               knitr_1.42                httpcode_0.3.0